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Abstract 

Native oyster populations in Chesapeake Bay have been the focus of three decades of restoration 
attempts, which have generally failed to rebuild the populations and oyster reef structure. Recent 
restoration successes and field experiments suggest that high-relief reefs offset heavy sedimentation 
and promote oyster survival, disease resistance and growth, in contrast to low-relief reefs which 
degrade in just a few years. These findings suggest the existence of alternative stable states in 
oyster reef populations. We developed a mathematical model consisting of three differential equa- 
tions that represent volumes of live oysters, dead oyster shells (= accreting reef), and sediment. 
Bifurcation analysis and numerical simulations demonstrated that multiple nonnegative equilibria 
can exist for live oyster, accreting reef and sediment volume at an ecologically reasonable range of 
parameter values; the initial height of oyster reefs determined which equilibrium was reached. This 
investigation thus provides a conceptual framework for alternative stable states in native oyster 
populations, and can be used as a tool to improve the likelihood of success in restoration efforts. 
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1. Introduction 



1.1. Decline and restoration of native oyster populations 

In the pas t century, the native Ch esapeake Bay oyster, Crassostrea virginica, a dominant ecosys- 
tem engineer (jCerco and Noel . 120071 ) . has dr opped to approximately 1% of its previous abundance 



due to overfishing and habitat degradation ( Rothschild et al. . 19941 ). Harvests peaked in 1884 at 



615,000 metric tons, but in 1992 the harvest was only 12,000 metric tons. Additionally, human 
activities on land increased the flow of sediment into the bay 's waters, which weakened phys i ologi- 
cal health, lowered fe cundity and raised mortality of oysters dNeweli Il988l : iRothschild et all Il994l : 



Lenihan et al 



199sh . Exacerbati ng the situation, the p hysical profile of reefs has been leveled 
by fishers exploiting oyster reefs ( Rothschild et al. . 19941 ). which places the oysters lower in the 
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water column wher e water flow is reduced and sedim ent accumulation rates are highest, thereby 
suffocating oysters ( Newell . 1988 : Lenihan et al. . 19991 ). 



Efforts to restore native oyster populations have been extensive but largely ineffectual (jRuesink et al. 



20051 ). However, a recent restoration effort in the Great Wicomico River has yielded promising re- 
sults. In 2004, the Army Corps of Engineers created nearly 40 ha of reef in the G reat Wicomico 



River (Chesapeake Bay) consisting of oyster shell planted at different reef heights (jSchulte et al 



20091 ). The field experiment featured high-relief reefs built at an average of 25 — 42 cm in height, 
low-relief reefs at 8 — 12 cm, and controls of unrestored bottom. After three years, in 2007, the 
higher reefs were considerably more successful. Mean oyster density was fourfold higher on high- 
relief reefs, about 1000 oysters per m 2 , than on low-relief reefs. Moreover, the high-relief reefs 
were robust in architecture and resistant to natural disturbances, whereas the low-relief reefs were 
heavily sedimented and apparently on a trajectory to a degraded state. These differing outcomes 
suggested the potential for alternative stable states driven by the initial condition of reef height. 



1.2. Alternative stable states 

The dramatic decline in the oyster population as well as the marked difference in success 
of the high-relief and low-relief reefs may be explained in the context of catastrophic shifts and 
bistability (i.e., alternative stable states). Some ecosystems ex hibit precipitous shif ts in state 



without correspondingly dramatic changes in external conditions dScheffer et all lioOlh. Gradual 



changes in the external conditions of a system produce correspondingly gradual alterations in state 
variables until sudden drastic changes occur in the state variable. After a degraded state has been 
reached, return to the external (i.e., environmental) conditions present before the change does 
not induce reversion of the system to the pre-shift state. These systems exhibit bistability, which 
is one man ifestation of the exis t ence of mul t iple equilibria for a range of constant environmental 
conditions dScheffer et all l200ll : bong et all 120021 : iGuilll . 120091 ) . 



Alternative (= multiple) stable states in communities have been tr iggered by environmental 
distu rbances in va rious ecosystems including coral re efs ( Hughes . 19941 ) , lakes (Carpenter et al. 



19991 ). grasslands ( Rietkerk and Van de Koppel . 1997 ). and kelp forests ( Konar and Estea . 20031 ) . 
In mollusks (e.g., mussels, clams, and oysters) alternative stable states have been documented in 



beds of the horse mussel Atrina zelandica in New Ze aland (Coco et "all. bOOfih and the blue mussel 



along the northeast Atlantic coast of North America ( Petraitis et al. . 20091 ). Despite the likelihood 



of multiple stable states in marine species, there have be en few mathematical models of this process, 
particularly for mollusks dPetraitis and Hoffmanl . lioiol). Moreover, for native oyster populations, 
the mathematical models of population dynamic s have usually emph asized linear interactions with- 



out the potential for alternative stable states ( Powell et al. . 20061 ). wher eas biologically realisti c 



models of oyster filtration have integrated nonlinear biological processes ( Cerco and Noel 2007 ). 
Consequently, our mathematical formulation represents an advance in the mathematical modeling 
of alternative stable states in the population dynamics of exploited marine species such as the 
oyster, and thereby advances the theoretical underpinnings for ecological restoration of marine 
species. 



1.3. Feedback mechanisms for alternative stable states 



A lternative st able states are generally due to one or more feedback mechanisms (jScheffer et al. 



2001 : Guill 20091 ). In the Chesapeake Bay, oysters filter phytoplankton (microscopic algae) and 
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sediment fl owing onto ree fs, which lowers turbidity levels and the incidence of low dissolved oxygen 
conditions (iNewell Il988l ). Historically, massive reductions in oyster biomass and degradation of 
the reef matrix contributed to increasing sediment in the water column and low dissolved oxygen on 
the bottom of the bay. Additionally, oysters encou nter a greater pr oportion of the sediment in the 
water column when they are closer to the bottom ( van Rijn . 19861 ). which occurs with reductions 
in the vertical relief of reefs. The sediment negatively affects oysters by causing them to expend 
energy to filter it, thereby increasing susceptibility to disease and mortality rates, while decreasing 
growth and reproduction. Raising the oysters in the water column can lead to h igher fecundity 



and decreased mortality fr om reductions in turbidity and elevated filtration rates (jLenihan et al. 



19991 : iRuesink et all . l2005h 



A well known ecosystem in which feedback mechanisms lead to alternative stable states is that 
of shallow lakes, which have ecosystem properties similar to those of Chesapeake Bay. In shal- 
low lakes, aquatic vegetation dampens resuspension of sediment, reduces nutrients in the water 
column, a nd pr ovides protection from fish predation for zooplankton that feed on phytoplankton 
(jSchefferl . l2009h . Fish control zooplankton and resuspend sediment and nutrients by disturbing 
the lake floor. However, aquatic vegetation and zooplankton have been depleted by herbicides and 
pesticides. The reduction in vegetation leads to an increase in nutrients in the water column which 
increases phytoplankton, which in turn feed on the nutrients. Declines of zooplankton result in 
unchecked growth of phytoplankton. A dramatic increase in phytoplankton precludes light from 
reaching the lake floor which causes the vegetation to decline further. This exposes zoopl ankton to 
incre ased predation, which in turn leads to lower predation pressure on phytoplankton (jSchefferl . 

To return the lakes to a state of high vegetation and controlled phytoplankton abundance, 
zooplankton populations must be rebuilt to a critical level such that they can control phytoplank- 
ton, and thus facilitate the profusion of vegetation that provides protection from fish predation. 
By temporarily reducing fish abundance, zooplankton can proliferate and feed on the abundant 
phytoplankton. Once phytoplankton are reduced, the vegetation can become re-established. The 
increased zooplankto n eventually al lows the fish population to rebound, restoring the system to the 
pristine stable state ( Scheffer . 20091 ). 



We hypothesize that the oyster reef system is somewhat analogous t o that of sha llow lakes. 
Oysters can control the volume of sediment while consuming phytoplankton ( Newell . 19881 ). but they 
must first be provided with op timal reef features, p articularly an elevated reef height. This enhances 



recruitment of young oysters ( Schulte et al. . 20091 ). which subsequently grow to a dense spawning 



stock that further filters phytoplankton and sediments from the water column and preclude the 
accumulation of sediments on the reef. Without the initial reef height it has been postulated that 
the low-relief reef structure inhibits high recruitment, resulting in sparsely distributed oysters that 
cannot keep pace with the sediments and phytoplankton, eventually leading to a heavily silted, 
degrading reef. We now provide the theoretical underpinnings for this interaction between oysters, 
reef height and sediment, and demonstrate that alternative stable states are feasible for oyster reef 
populations. 



1.4- Objectives 

We construct an ordinary differential equation model to demonstrate that the mechanisms 
involved in the interaction of oysters, reef height and sediment produce bist a bility, which provides 



an explanation for the success of high-relief reefs and failure of low-relief reefs( Schulte et al. . 20091 ) 
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in the context of multiple stable states. Bifurcation theory is used to identify parameter ranges 
that produce bistability in the model. The bistable structure is sensitive to the initial values 
of the system. A small perturbation of the initial value can change the eventual outcome from 
one stable state to a different one, which is termed a "hair-trigger effect". In gener al the basins of 
attra ction of the two stable states are only separated by a surface in the phase space (jJiang and Sh.il . 
We demonstrate our model's sensitivity to initial conditions via numerical simulations. The 



mathematical model of oyster and sediment is presented and motivated in Section [2j Analytic 
bifurcation results are given in Section \3.1\ and numerical simulations for specific parameter values 
are in Section 13.21 In Section U] we discuss the implications for alternative stable states in oysters 
for the likelihood of success of native oyster restoration, with emphasis on the Chesapeake Bay 
ecosystem. 



2. Methods: mathematical model 

2.1. Basic elements of the model 

We model the rate of change of live oysters, dead oyster shells, and deposited sediment volume 
with respect to time, t, measured in years. Our state variables are measured in volume per m 2 of 
sea floor, so they can be converted to heights measured from the sea floor by dividing by a unit 
area. The feedback interaction between live oysters and sediment occurs as follows. Live oysters 
follow logistic growth but are negatively affected by sediment volume. We introduce a function to 
represent the proportion of oysters above the level of suffocating sediment. The change in dead 
oyster shell volume is due to the death of live oysters minus a degradation rate proportional to 
dead oyster volume. The volume of sediment deposited on the reef depends on the position of the 
reef in the water column and on filtration by live oysters. We next give a detailed derivation for 
each differential equation. Variables are summarized in Table Q] and parameters in Table [2j 



Variable 


Description 


Units 


t 


time 


years 


0(t) 


live oysters 


q 

m° 


B(t) 


dead oyster shells 


•i 

m° 


S(t) 


sediment on reef 


m i 



Table 1: Model variables. 



2.2. Proportion of oysters unaffected by sediment 

We introduce a function / that represents the proportion of oysters not affected by sediment. 
The input, do, represents the volume of the live oysters (O) and dead oyster shells (B) not affected 
by sediment (S). Hence we define 

(1) do = O + B - S. 

Note that the volume is essentially height multiplied by a unit area, so do represents the height 
of the oysters above the sediment. In this we make a simplifying assumption that there is a 1:1 
relation between oyster or shell volume and sediment volume. Future work will model in more 
detail how sediment volume is affected by reef geometry. 
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The proportion of oysters not affected by sediment is an increasing function of do € (—00,00) 
with a sigmoid shape bounded by and 1. We denote this function by /i(do). When do = O, the 
dead oyster shells are covered by sediment but the live oysters are not, i.e., when B = S. In this 
situation, live oysters are not affected by sediment and fi(0) ~ 1. When do = 0, both live and 
dead oysters are covered by sediment, and /i(0) ~ 0. As do approaches the positive and negative 
limits, the approximations are equalities. We can assume that /i(do) is a sigmoid function with the 
only reflection point at the midpoint of the interval [0, O] and /i(0/2) = 1/2. By the replacement 
d = do — 0/2, the function /i(do) becomes a function of: 

(2) d = do -0 = ° +B -S. 
Then we define f(d) = f±(do — 0/2). In summary, f = f(d) satisfies 

(3) /'(d) > 0, /(0) = \, lim f(d) = 0, and lim /(d) = 1. 

Z a— ¥— 00 a— ¥00 

The function / was devised to make /(d) ~ when d < —k\/2 and /(d) ~ 1 when d > fci/2, where 
k\ is a parameter proportional to the carrying capacity of the live oysters. Thus, 0/2 oysters must 
be covered by sediment before their performance is affected by sediment; i.e., we assume that more 
than half of any individual oyster must be covered by sediment for its performance to suffer. 



2.3. Live oyster volume 

The change in live oyster volume 0(t) is represented by the differential equation 

(4) ^ = rOf(d) (l " f ) - »f(d)0 - 6(1 - f(d))0. 

Here the live oysters are assumed to follow logistic growth, K is the carrying capacity, and r 
represents the intrinsic rate of increase. The oyster population increases at a negative density- 
dependent rate until the population reaches K. In the second term of the equation, fi represents 
mortality due to predation and disease. Both terms are scaled by /(d) because oysters covered in 
sediment do not reproduce, and are assumed to die due to suffocation by sediment rather than by 
predation or disease. The third term represents the decrease in live oyster volume as a result of 
sediment; e is the death rate of oysters covered by sediment. As /(d) goes to 1, this term goes to 
zero. However, as /(d) becomes smaller, the term begins to exert a greater effect. 



2-4- Dead oyster (shell) volume 

The second differential equation represents the change of dead oyster shell volume B(t) as 

(5) ^ =fim0 + e(l-f(d))0- 7 B. 

The first two terms are directly from the death of live oysters in Equation , and the third term is 
the loss of dead oyster volume due to degradation of shell. This loss is proportional to the volume 
of dead oysters at the rate 7. Note that the term r0 2 f(d)/K in is not included in ([5]) as it is a 
reduction in growth rate rather than a loss term, and it does not increase the dead oyster volume 
as do the other two terms. 
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2.5. Sediment accumulation 



The system of differential equations is completed by a third equation describing the change of 
sediment volume S(t) as: 



(6) 



dS_ 
~dt 



-PS + Cge 



FO 
' °9 . 



Here the first term is the rate of sediment erosion, which is proportional to the volume of deposited 
sediment with a rate /3; the second term is the rate of sediment deposition. The sediment deposition 
rate in the absence of oysters is Cg, where C is a maximum possible deposition rate and g is a 
modification that depends on reef height O + B. The deposition rate i s at its maxim um at the sea 
floor, and decreases as the reef height in the water column increases riNielseilfliii ). Hence with 
the reef height represented by x = O + B, we assume that the function g(x) = g(0 + B) satisfies 



(7) 



g(0) = 1, g'{x) < for x > 0, and lim g(x) = 0. 



In this mo del formulation we assu me that biodeposition is a constant, minor fraction of sediment 
deposition ( Cerco and Noel . 20071 ). and which would not alter the qualitative results of the mod- 
eling. Other processes that are not simulated explicitly, such as increased deposition due to feces 
and pseudofeces at the bottom, are parameterized by the erosion and burial rates. In future, more 
complex formulations we will integrate biodeposition into the model to generate results that are 
more accurate quantitatively. 

In the presence of live oysters, the deposition term should be reduced by a multiplicative factor 
due to filtration. The filtration rate F per unit oyster volume depends on the height-dependent 
sediment concentration, which is proportional to Cg. The rate F scales linearly with Cg when Cg 
is small, reaches a peak Fq at some optimal sediment concentration, and beyond this threshold, it 
decreases as oyster gills become increasingly clogged (jJordanl . Il987l ). Hence F = F(y) = F(Cg) 
satisfies 



(8) 



-F(O) = 0; lim F(y) = 0; and there exists yo > such that 

y-¥co 

F'(y) > for < y < y , F'(y) < for y > y , and F(y ) 



Both g(x) and F{y) are positive and continuously differentiable functions on [0, oo). 

The deriva tion of (|6l) beg ins with a mass balance of sediment deposition in a control volume 
at the bottom (jChapral . Il997t IH l2008h : 



(9) 



d(VC s 
dt 



Q s AC b - V e AC s - V b AC s - V f O, 



where £l s is settling velocity (m/day), V e is erosion velocity (m/day), V b is burial velocity (m/day) of 
sediment that is consolidated in the bottom and lost from the system, C s is sediment concentration 
at the bottom (g/m s ), C b is the sediment concentration above the bottom sediment layer (5/m 3 ), 
V is the control volume (m 3 ), A is the surface area (m 2 ) of the control volume, Vf is the filtration 
rate by oysters (g/(m 3 ■ day)), and O is live oyster volume (m 3 ). 

We introduce a mean erosion rate v e (1/day) and burial rate v b (1/day), which can be considered 
to be a re-scaling of the erosion and burial velocities V e and V b by the mean depth as follows: 

d{VC s 



(10) 



dt 



n s AC b - v e VC s - v b VC s - VfO. 
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The bottom se diment concen tration can be expressed as sediment density p and the porosity 
as C s = p(l — 4>) ( Chapral . 19971 ). We assume that the mean porosity, 4>, is a constant and divide 



V(l 

(|10p by p(l — 4>). Let S = = — be the sediment volume with porosity normalized by the 

1- (J) 

mean porosity to quantify the volume of unconsolidated sediment deposition at the bottom, and 

tt s A V f . . 

let u> s = — =- and vr = — — — =r. Then (|10j) can be written as 

p(l - <p) p(l - <f>) 

(11) ^- =u s C b -v e S -v b S -VfO. 



Rearrangement of the terms in Equation (jlip yields 
(12) | = „, Ct ( 1 _^|)_ ( „ e + „ t ) S . 

Here the first term represents the sediment deposition modified by oyster filtration, and the second 
term is the combined loss of sediment due to erosion or burial. Because the modified deposition 
term could be negative for large values of O, we instead replace the decreasing linear function 

VfO . f—VfO\ vrO 
1 — — by a decreasing nonlinear function exp — —— whose linearization is 1 —pr- Now 

uj s C b \ UJ S C b J LO S C b 

we rename 

cosC b = Cg, and vj = F, 

where C is a constant representing the maximum deposition rate, and g is a decreasing function of 
O + B with maximum g(0) = 1; then we obtain ([6]). The estimates of C and g come from data of 
ujs an d C b ; information on Vf can be used to determine the form and parameter values of F; and, 
v e and v b determine (5. 

Equation ([6]) has properties reflecting the desired qualities of the system. When O — ¥ 0, 
S ~ Cg — /3 S (deposition and erosion without oysters). When O is small, from a Taylor expansion 
of the exponential we have S ~ Cg — FO — f3S so the deposition rate is reduced linearly by the 
amount FO that oysters filter out of the system. When O — > oo, S ~ — fUS (no deposition, only 
erosion). Finally, when the filtration rate per unit oyster volume F — > 0, either because there is 
too little sediment to filter or the oysters are being choked by it, S ~ Cg — j3S. 



2.6. Full model 

In summary, we propose the following set of differential equations to model oyster population, 
oyster reef, and sediment volumes: 

(13) ^ = rOf(d) (l - x) " Vf(d)0 - 6(1 - f(d))0, 

(14) d ^ = llf{d)0 + e(l-f{d))0- 1 B, 

dS FO 

(15) _ = _ /3< s + C 5 e-^, 

where the quantities d, f(d), g = g(0 + B) and F = F(Cg) satisfy ([3]), ([7]) and (jSJ), respectively. A 
set of functions satisfying these conditions will be given in Section 4. The parameters of the system 
(|13p -(|15p are summarized in Table [2j The last two parameters h and r\ are specific to the choice of 
/ and g, which will be explained in Section [3l 



7 



Parameter 


Meaning 


Units 


Value 


Reference 


r 


instantaneous rate of increase 


year -1 


0.7-1.3 


1 


K 


oyster carrying capacity 


m i 


0.1 - 0.3 


2 


/i 


mortality rate due to 
predation and disease 


year -1 


0.2 - 0.6 


3 


e 


mortality rate due to sediment 


year -1 


0.94 


4 


7 


oyster shell degradation rate 


year -1 


0.5 - 0.9 


4 


^0 


maximum sediment filtration rate 


year -1 


1 


5 


C 


maximum sediment deposition rate 


m 6 ■ year -1 


0.04 - 0.08 


6, 7 


yo 


sediment amount where 
the filtration is maximum 


year ■ m -i 


0.02 


5 




sediment erosion rate 


m -i 


0.02 - 0.04 


6,7 


h 


scaling factor 


m 


10 - 30 




V 


decay rate of sediment deposition 
on the reef height 




3.33 


8 



T able 2: Table of para met ers in the equation s. R eferences: 1 |U.S. Army Corps of Eng ineers Norfolk District , 2009), 
2 (ISchulte et all l2009h, 3 dPowell et al.l , l2009l ), 4 (|Smith et al.l . l2005l) . 5 (|Jordanlll987l ), 6 (|Kniskern and Kuehlll2003h , 
7 (|chapral,ll997T l, 8 (|van Riinl , ll986l ). 



3. Results 



3.1. Bifurcation analysis 

We are interested in determining the values of parameters and state variables at which the 
change in the state variables is equal to zero (i.e. a steady state reef-sediment system). We 
consider the equilibrium solutions of our model f)13[) -(|15 [) . which satisfy 



(16) 

(17) 
(18) 



= rOf(d) 1 



O 
K 



M /(d)0 - e(l - f(d))0, 



= Hf(d)O + e(l-f(d))O-iB, 

FO 

= -pS + Cge~^ . 



A trivial solution of f)16|) — (|18p where O = B = and S = C / (3 is an equilibrium solution 
representing the extinction of the oyster population and the accumulation of sediment limited only 
by erosion. We will now solve the system in search of nontrivial solutions where O > 0, B > and 
S > 0. 



We describe a procedure of reducing the equations (|T6j) - (|T8|) to a single equation. From (fT6|) . 
we obtain (assuming O > 0) 



(19) f(d) 
and similarly from (|17p . we obtain 

(20) f{d) 



eK 



K(r — fi + e) — rO ' 
^B-eO 
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From (|19p and (|20p . we can have an equation of O and B only: 

(21) * = 

V ' K(r - fi + e) -rO (/x-e)0' 

and B can be solved from (1211) as 



(22) 5 - 7 [^(r-^ + e )-rO] =B(0) - 

One can solve 5 from (1181): 



(23) S = 5(0,S) = ^e- F ' / cr », 

where 5 depends on O + B and i* 1 depends on g. Now the substitution of (|22p and (|23p into ([19 
results in an implicit equation of O only: 

(24) f (° + B(0)-S(0,B(0))^ 



2 v y s ' v ") Kir - /x + e) -rO 



Hence for a fixed set of parameters, any root O* > of (|24p corresponds to an equilibrium point 
(0*,B(0*),S(0*,B(0*))) of (fT6j) - ([T8j) . While direct analysis of {24]) is not simple due to the 
complicated definitions of 0(B) and S(0,B), numerical calculation of (|24p is relatively easy with 
a symbolic mathematics software. 

We define the functions on the left and right hand side of (124)) to be 
(25) L(0) = f ( £ + fl(O) - S(0, B(0))) , 



(26) R(0) 



2 



K(r — jx + e) — rO 



From (I24p . intersection points of the graphs of -L(O) and R(0) are equilibrium points. We observe 
that L(0) is bounded by 1 and R(0) is unbounded as O — > = K[r — fx + e)/r. Thus if 
L(0) > -R(O), then (|24h has at least one root with positive O from the intermediate- value theorem; 
and if L(0) < R(0), then (|24f) may have no or two zeros for most cases, or one zero in the case that 
L(0) and R(0) are tangent to each other. From (|16p we see that 

L(0) = /(-C//3), and i?(0) = 



- \x + e 

From a different point of view, one can consider the equations (|16p - (|18p with a bifurcation 
analysis and linearization. Linearizing (|16p - (fl~8l) at the trivial equilibrium (0,B,S) = (0, 0, C//3), 
we obtain the Jacobian matrix to be 

/ f(-C/P)(r-n + e)-e 
(27) J(0,0,C//?)= /(_cr//9)( M -e) + e -7 

V CV(0)-F(C) C 5 '(0) -/3 

Since the Jacobian matrix J(0, 0, C//3) is lower triangular, the three diagonal entries are eigenvalues. 
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We discuss the equilibrium problem in several cases: 

Case 1 : If r < fx (the birth rate smaller than the natural death rate), then the trivial one (0, 0, C / (5) 
is the only equilibrium. In this case O' < 0, thus the live oysters are destined to go extinct. So the 
equilibrium (0,B,S) = (0,0, C/ (3) is globally asymptotically stable. 

Case 2 : If \x < r < /i + e (the birth rate larger than the natural death rate, but smaller than the 
combined death rate due to natural cause and due to sediment), positive equilibrium points are 
possible. We notice that f(0) = 1/2 from ([3]), so r < fi + e is equivalent to 

f(0)(r-n + e)<e. 

This implies that for any C > 0, /(— C//3)(r — /i + e) — e < /(0)(r — fi + e) — e<0 since / is an 
increasing function. Thus the trivial equilibrium (0, 0, C//3) is locally stable for any C > 0. 

For a critical value r* G + e), if r > r*, then equation (f2i|) has exactly two positive roots 
when C = if f(d) is a concave function for d > 0. And for a fixed r value satisfying r* < r < /x + e, 
there exists another critical value C*(r) > such that equation (|24p has exactly two positive roots 
for < C < C*(r). When the maximum sediment deposition rate C is large, the system (jl3l) -(jl5[) 
can only have the trivial equilibrium. Hence the parameter region for existence of two positive 
equilibria when /i<r</i + eis {(r, C) : r* < r < /j, + e,0 < C < C*(r)}, given that all other 
parameters are fixed. 

Case 3 : If r > fi + e (the birth rate larger than the combined death rate due to predation, disease 
and sediment), then 

f(0)(r-n + e)>€. 

The monotonicity of / assumed in (|3|) implies that there exists a unique C*(r) > such that 
f(-C/f3)(r - n + e) - e > for C > C*{r), and f(-C/f3)(r - fj, + e) - e < for C < C*{r). Thus 
in this case, the equilibrium (0, 0,C//3) is locally stable when C > C*(r), and it is unstable when 
0<C<C*(r). 

The critical value C = C*(r) is a bifurcation point where a branch of nontrivial equilibrium 
points emanates from the line of trivial equilibria (C,0,B,S) = (C, 0, 0, C/(3). The bifurcation is 
backward if the bifurcating equilibria are unstable, otherwise it is forward: the bifurcating equilibria 
are stable (Fig. [1]). When the bifurcation is forward, a unique equilibrium exists for < C < C*(r), 
and there is no nontrivial equilibrium for C > C*(r). Conversely, if the bifurcation is backward, 
then there is a range of values of C > C*(r) for which the system has two nontrivial equilibria. If 
the sediment volume is too large, there is a largest value C = C*(r) such that the system only has 
the trivial equilibrium when C > C*(r). The parameter region for two positive equilibria occurs 
when r > /i + e is {(r, C) : r > /i + e, C*(r) < C < C*(r)}, given that all other parameters are fixed. 
Note that the bifurcation occurring at C = C^r) is a transcritical one. The terms "backward" 
and "forward" are abo ut the branch of positive equilibria, an d they are similar to the ones used in 
epidemic models (see ( Hadeler and Van den Driessche . 19971 ).) 



In the Appendix, we show that the direction of the branch of bifurcating positive equilibria is 
determined by 

(28) i = Ml - e + ^ K _ lt%m + M _ 4 (cm - no + c ' mcr 



2r n-C/P) 7 7/3 V ' 7(r-/« + t) 

If / < then the bifurcation is forward; if / > then it is backward which implies bistable 
parameter ranges. 
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Figure 1: Schematic Diagrams of the Forward and Backward Bifurcations. The horizontal direction is the parameter, 
and the vertical direction is the state variable. Stable equilibria: thick curves; unstable equilibria: thin curves, (left) 
forward bifurcation; (right) backward bifurcation. 

Therefore the question of bistability for this case is reduced to whether I > 0. Notice that 
c/'(0) < and F(cJ) > 0, hence 

= A(r - , + XjK _ ,K / _ N > 

2r 7 7P V 7( r _ M + e ) / 

and the positivity of / = Ji — 12 depends on the competition between I\ and I2 = -jf, — * , > 0. 
Here I2 only depends on the form of / and C*//3, while I\ depends on many other parameters. 

T — LJL ~\~ 6 

Notice that C* is determined by / and . If we fix the values of r,fi,e, and ft, then one 

e 

can increase I\ to generate bistability by (i) increasing carrying capacity K; (ii) decreasing the 
oyster shell degradation rate 7; (iii) increasing |</(0)|, the decay rate of sediment deposition with 
reef height; or (iv) increasing F(C*), which represents oyster filtration efficiency. 

3.2. Numerical simulations 

To illustrate our results numerically, we define the functions f{d), g(x), and F(y), which satisfy 
the desired conditions ([3]), ((7J), (jSJ): 

fW = TTT^> d=^ + B-s, 

(29) g(x) = e-f, (g(0 + B) = e -^°+ B )), 

F(y) = ^ye^ ~y^y°, (F(Cg) = ^£l e (yo-Cg)/yo^ 

yo yo 

Here h > is a parameter that adjusts the shape of the sigmoid function /. For larger h, the 
function / has a narrower transition where the function value jumps from to 1. In the definition 
of g, r\ is the decay rate of the exponential function. The per volume filtration rate F = F{Cg) 
is a function of Cg and is of the Ricker type, where Fq represents the maximum filtration rate, 
achieved at y = yo (Fig. [5]). 

With the nonlinear functions /, g and F defined as in (|29p . we assume the following set of 
parameters: 
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Cg 

Figure 2: Graphs of functions /(d), g{x) and F(y) defined in (f29)l . (left) /(d) with h — 20; (middle) g(x) with 
i) = 10/3; (right) F(y) with y = 0.02 and F = 1. 



Parameter 


r 


if 


M 


e 


7 


»7 


yo 


5 





h 


C 


Value 


1 


0.3 


0.4 


0.94 


0.7 


3.33 


0.02 


1 


0.01 


20 


0.02 



Table 3: A sample set of reasonable parameters 



For the parameters given by Tabled there are two positive equilibrium points (0i,5i,Si) = 
(0.0566,0.0456,0.0326) and (0 2 ,B 2 ,S 2 ) = (0.1736,0.1022,1.0645 x 10~ 7 ); thus, the parameter set 
in Table[3]is in the bistable region. Freeing the parameter C gives a bifurcation diagram (see Fig. [3]) 
with two positive equilibria for all < C < C*, where C* ~ 0.078 is a saddle-node bifurcation 
point where the curve bends back. This bifurcation diagram confirms the description we give in 
Section [3] Case 2, as 0.4 = /x<r = l</i + e = 1.34. 

We use numerical simulation to examine the bistable dynamics of (|16jl - (jl8j l with parameters 
given in Table [3l We use the initial value of O(0) = 0.01 and S(0) = 0.01; i.e., there is a small 
amount of live oyster and also a small amount of sediment initially. We chose several different values 
of 5(0): 5(0) = 0.20, 0.10, 0.12 and 0.11 (Fig. gj). For larger 5(0), the oyster population survives 
and reaches the large stable equilibrium (Oi, Bi, Si), whereas the smaller 5(0) will drive the oyster 
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Figure 3: Bifurcation diagram of (|16|) - (I18| I. with parameters given in Table [3] Here the horizontal axis is sediment 
deposition rate C, and the vertical axis is (upper left) live oyster volume O, (upper right) dead shell volume B and 
(lower) deposited sediment S. The point labeled 'LP' in the interior is the saddle-node bifurcation point, and the 
point labeled 'H' is a neutral saddle point. The bifurcation diagram is generated by Matlab package MatCont. 

population to local extinction. The critical level of initial reef height B(0) is between 0.11 and 0.12. 
The transient dynamics with -B(O) = 0.12 and B(0) = 0.11 when < t < 10 (Fig. [5]) demonstrates 
that at the higher reef height, live oyster volume increases and eventually curbs sediment volume 
to a very small level due to oyster filtration and reef height. The slightly lower reef does not permit 
live oysters to filter the sediment sufficiently, and the sediment eventually covers both live and 
dead oysters. The specific reef heights (i.e., 0.11 and 0.12) that discriminate the two trajectories 
towards stability are likely to differ depending on natural variation in other parameters of the three 
equations, such that they should not be viewed as rigid values under all conditions. Rather, the 
key point is that a slight shift in initial conditions can drive the system towards two dramatically 
different trajectories. 

The bistability displayed with parameter values given in Table [3] is not anomalous (Section 
[3|); the parameter range for bistability is robust. For example, given all other parameter values 
as in Table [3] except the oyster growth rate r and sediment deposition rate C, then the range of 
parameters (r, C) to produce bistability has been shown to be 

{(r,C) : r* < r < /i + e,0 < C < C*(r)} and {(r, C) : r > fi + e, C*(r) < C < C*(r)}. 

The graphs of C*(r) (Fig. [U upper curve) and C*(r) (Fig. [HJ lower curve) are monotonically 
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Figure 4: Numerical solution (0 < t < 50) of (I16I) -(I18 | ) with parameters given in Table [3] for various initial reef 
heights. For all cases O(0) = 0.01 and 5(0) = 0.01. (upper left) B(0) = 0.20; (upper right) 5(0) = 0.10; (lower left) 
B(0) = 0.12; (lower right) B(0) = 0.11. 



increasing functions of r. The first critical value r = r* ~ 0.739 is where the curve C = C*(r) 
emerges from C = 0, and the second critical value r = \i + e = 1.34. The region between the 
two curves is the "bistable regime", and the chosen parameter value (r, C) = (1,0.02) is in that 
region. The region above the bistable one is the "extinction regime" where the trivial equilibrium 
(0, 0, C//3) is globally asymptotically stable; the region below the bistable one (which is very small 
and can only be detected on the log-plot) is the "persistent regime" where the oysters will persist 
irrespective of initial reef height. Note that the estimated growth rate range 0.7 < r < 1.3 is 
dominated by the bistable regime. 

The system can shift from Case 2 to Case 3, where a transcritical bifurcation occurs for a 
positive C* (r) (Fig. [7|) when r is large or when e is low. 
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Figure 5: Numerical solution (0 < t < 10) of (I16I) - (I18[) with parameters given in Table [3] for various initial reef 
heights, (left) Initial value O(0) = 0.01, 5(0) = 0.12 and 5(0) = 0.01; (right) Initial value O(0) = 0.01, B(0) = 0.11 
and 5(0) = 0.01 



4. Conclusions 

4-1. Key findings 

We constructed an ordinary differential equation model for the volumes of live oysters, shell 
from dead oysters, and accumulated sediment on an oyster reef. Feedback interactions between 
the oysters and the sediment occurred such that low-relief (i.e., low vertical height) reefs were 
eventually choked by sediment, whereas high-relief reefs had less sediment deposition due to their 
height off the sea floor and filtration of the sediment by live oysters. Using sediment deposition rate 
as a bifurcation parameter, we have shown that the oyster-free (i.e., degraded reef) equilibrium 
point is stable for sufficiently high sediment influx. In contrast, it is unstable at lower sediment 
deposition rates, and in that case the system has at least one equilibrium point with live oysters 
and a persistent reef matrix. Furthermore, the system has a backward bifurcation in sediment 
concentration when certain criteria are met, whereby there are two stable states, one with live 
oysters and persistent reef, and one without live oysters or shell reef. Bistability becomes more 
likely when the oyster reefs grow higher from the sea floor and degrade less quickly, and when 
the reefs more strongly reduce sediment deposition either through reef height or through enhanced 
filtering of sediment by live oysters. 

We also observed bistability in numerical simulations for physically realistic parameter values. 
In this case, the long-term fate of a reef depended on its initial conditions. In particular, an initial 
low volume of dead shells (i.e., reef height) led to extinction of the live oysters and degradation 
of the reef, whereas an initial high volume of dead shells led to a stable steady state with living 
oysters. These res ults are analogous to the experimental field results for low-relief and high-relief 
reefs, respectively ( Schulte et al. . 20091 ). This study is therefore the first to provide a theoretical 



foundation for the existence of bistability in oyster reefs. 
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Figure 6: The graph of C*(r) (curve of saddle-node bifurcation points) and C*(r) (curve of transcritical bifurcation 
points), with parameters given in Table [3] except r and C. Here the horizontal axis is oyster growth rate r, and the 
vertical axis is (left) log(C), (right) C. 



4-2. Model assumptions and caveats 



The assumptions made in our model are the following. Oysters grow logistically and die linearly 
when above the sediment. They die at a constant rate per oyster volume when buried by sediment. 
To quantify the extent to which oysters are above the sediment, we use a monotonically increasing 
function / € [0, 1] that depends on the difference between the reef height and the sediment height. 
Dead shells in the reef are degraded at a constant rate per dead reef volume. Sediment is eroded 
at a constant rate per sediment volume. It is deposited at a rate that decreases with the height 
of the reef in the water column and decreases exponentially with the rate at which oysters filter 
sediment. The filtration is assumed to be proportional to the oyster volume and to a filtration 
rate that has a single maximum at some optimal sediment concentration. These assumptions 
are sufficient to generate regions of bistability. For the numerical studies, we further assumed a 
particular sigmoidal form for /, we assumed that sediment concentration decreases exponentially 
with height in the water column, and we assumed that the filtration rate is a Ricker-type function 
of the sediment concentration. 

For direct co mparison of our res ults with experimental field data from natural systems (e.g., 
Chesapeake Bay ( Schulte et al. . 20091 )). reliable estimates are needed for oyster growth and param- 
eters related to the interaction between oysters and sediment. Further refinements to the model 
may include seasonal effects, spatial effects such as sources of oyster larvae from distant reefs, and a 
more accurate representation of reef geometry, including how sediment intercalates between shells. 



4-3. Relevance for oyster reef restoration 

In previous mathematical and c onceptual models of oyster reef dynamics, the possibility of 
alternative stable states was ignored (jPowell et all 120061: iMann and Powelll . 120071 ). despite the evi- 
dence for nonlinear processes in oyster ecology (jCerco and Noell . 120071 ) . This omission promoted a 
narrow focus regarding the optimal reef architecture in oyster restoration efforts within Chesapeake 
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Bay, despite the highly variable environmental conditions (e.g., sediment deposition) throughout 
the bay's wate r s. Th e collective findings of our theoretical analysis and previous field experiments 
(|Schulte et all 120091 ) demonstrate that the optimal oyster reef architecture will differ based on 
environmental conditions, and that initial reef height is a critical feature of reef architecture, which 
determines the eventual persistence or degradation of constructed oyster reefs. These results are 
not surprising given the existen ce of alternative s table states in othe r bivalve mollu s ks su ch as the 
horse mussel Atrina zelandica ( Coco et al. . 2006 ) and blue mussel ( Petraitis et al. . 20091 ). and in 
a diverse suite of ecosystems ( Scheffer et al. . 200ll ). and underscore the need to consider the phe- 
nomenon of alternative stable states in ecological restoration. This investigation thus provides a 
conceptual framework for alternative stable states in native oyster populations, and can be used as 
a tool to improve the likelihood of success in ecological restoration. 
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Appendix: Calculation of the linearization and stability 



In this section we explain the method of bifurcation of equilibrium points from a known branch 
of trivial equilibria. It is well-kn own as "bifurcation from a simple eigenvalue" in the studies of 



analy tical bifurcation theory (see (jCrandall and Rabinowita . 11971 



Jiang and Shi . 20091 ; Liu et al 



20071 )). Here we apply this powerful method to a finite-dimensional problem. 

Consider a smooth mapping F = F(X,u) : R X U — > W 1 where U is an open subset of W 1 , 
n > 1, A is a parameter and u is the state variable. We consider the equilibrium problem 

(30) F(A,«) = 0. 

Assume that a trivial solution is known. That is, there exists uq € U so that F(X, uq) = for all 
AgR. So {(A, no) : A E R} is a line of trivial solutions of (f30|) . 

The linearization of F with respect to u is represented by the Jacobian matrix: F u = (Jy = 
djFi), where Jij is the entry of F u at row % and column j, and 



dF(\,u) 
d Uj 



1 < i,j < n, 



is the partial derivative. Note that F = (F±,F2, ■ ■ ■ ,F n ) and u = (u\,U2, • • • ,u n ) are both vectors 
in W 1 . Similarly the second derivative of F on u is expressed as a 3-dimensional matroid F uu = 
(K ijk = d jk Fi), where 

9 2 F(A,n) 



dujduk ' 



1 < k < n, 



is the second order partial derivative. Also the mixed derivative F\ u 

d 2 F i {\,u) 



d\jFi) where 



dxjF 



dujdX 



1 < i, j < n. 



We notice that F u defines a linear operator R n — > M. n with matrix multiplication and so does F\ u . 
F uu defines a bilinear operator IR n x R n — >• ]R n which can be expressed as 



Fuu [(-^1 5 ' ' ' j %n) 



(Vl, ■ ■ ■ , Vn)] = K ljkXjVk, 
3,k 



^ ] K-njk%jyk)- 



Finally for a linear operator L : M n — > M n , we use N(L) and R(L) to denote the null space and the 
range of L; and we use (x,y) to denote the standard dot product of x,y 6 M n . 



Now we are ready to state a bifurcation theorem due to Crandall and Rabinowitz (jCrandall and Rabinowita . 

Il97lh (here we only state a special case): 

Theorem 4.1. Let F :Exf/ — >■ M n 6e tozce continuously differentiable, where U is an open subset 
o/M n . Suppose that F(\,uo) = /or AgK, and erf (Ao,«o)> F satisfies 

(Fl) dimN(F u (\Q,uo)) = codimR(F u (\o,uo)) = 1, and N(F u (Xq, uq)) = Span{wo} ; 
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(F3) F Xu (X ,u )[w ] £ R(F u (X ,u )). 



Then the solutions of (|30j) near (Ao,no) consists precisely of the curves u = uq and (X(s),u(s)), 
s €. I = (—5,5), where (X(s),u(s)) are continuously differ entiable functions such that A(0) = Ao, 
u(0) = uq, u'(0) = wq. Moreover 

(l,F uu (Xo,uo)[wo,w ]) 



(31) A'(0) 
where I € R n satisfying R(F u (Xq,uq)) 



2{l,Fxu(Xo,u )[wo\) 
{yeR n :(l,y)=0}. 



In simpler terms, at a bifurcation point A = Ao, the Jabobian F u has zero as an eigenvalue; 
(Fl) means that zero is a simple eigenvalue of F u , which means that the eigenspace of F u is one- 
dimensional, and the range of F u is (n — l)-dimensional (called, codimension one); (F3) means 
that F\ u [wo] does not belong to the range of F u , where wq is any nonzero eigenvector. Once 
these conditions are satisfied, then there is a curve of solutions bifurcating from the branch of 
trivial solutions. The formula of A'(0) is useful for determining the direction of the bifurcation 
(forward /backward) . 

To apply the above abstract theory to our problem, we notice that the trivial equilibrium 
(C, O, B, S) = (C, 0, 0, C / f3) is not constant for C. Hence we make a change of variable z = S — C/0 
then (O, B, z) = (0, 0, 0) is a constant solution. We define 



(32) 



G(C,0,B,z) 



I 



\ 



rOf(d) 



el 



rf(d)— + f if(d)0 

FO 
' Cg 



fxf(d)0 
jB + e(l 



~ f(d))0 
f(d))0 



c 



Cge c n - (5z 

where d = X(0/2 + B) — z — C/(5, and definitions of /, k, g, F are same as before. Let u = (O, B, z). 
Then G U (C, 0,0, 0) is the same as ([27|) . At C = C*, G U (C*, 0,0,0) can be written as 



(33) 



L = G U (C„ 0,0,0) 





cr 



We take the eigenvector of L to be wo 

WQ2 



r - /i + e 
V C*(/(0) - F(C„) 

= (1,^02,^03) where 
cr 





-7 
C*g'(0) 



\ 



WQ3 



7(r-^ + e)' 
C*g'(0) - F(C* 



+ 



C*g'(0)er 



(3<y(r - n + e) ' 

one can see that the range of L is {(0, y, z) G M 3 } which is two-dimensional, so we can take the 
vector I to be (1, 0, 0). A vector v does not belong to the range of L if the first entry of v is not zero. 
So to apply (|3ip . we only need to calculate the derivatives from the first equation of the system. 

We can calculate that (/, G^ U (C*, 0, 0, 0)[wo]) = — /'(— C*//3)(r — fx + ej/P < 0, and with a more 
tedious calculation, we find that 

(Z,G U u(*,0,0,0)[u;o, w ]) 



2r 

T 



f'(-c*/p) 



X(r-n + e)k f(-C*/0) Xek ek 



2r 



+ 



f'(-CJP) 7 7/3 



C*g'(0) - F(C») + 



CW(0)er • 
7(r-// + e), 
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Hence combining all the calculations, we obtain the direction of the bifurcating curve as 



C'(0) 



(J,G uu (C«,0,0,0)[wo,w ]) 
2(l,Gxu(C*,0,0,0)[w }} 



k(r - // + e) ' 



rl 



References 

References 

Carpenter, S., Ludwig, D., Brock, W., 1999. Management of eutrophication for lakes subject to 
potentially irreversible change. Ecological Applications 9 (3), 751-771. 

Cerco, C, Noel, M., 2007. Can oyster restoration reverse cultural eutrophication in Chesapeake 
Bay? Estuaries and Coasts 30 (2), 331-343. 

Chapra, S., 1997. Surface Water-Quality Modeling. McGraw-Hill. 

Coco, C, Thrush, S., Green, M., Hewitt, J., 2006. Feedbacks between bivalve density, flow, and 
suspended sediment concentration on patch stable states. Ecology 87 (11), 2862-2870. 

Crandall, M. C, Rabinowitz, P. H., 1971. Bifurcation from simple eigenvalues. Jour. Func. Anal. 
8, 321-340. 

Dong, Q., McCormick, P. V., Sklarb, F. H., DeAngelis, D. L., 2002. Structural instability, multiple 
stable states, and hysteresis in periphyton driven by phosphorus enrichment in the everglades. 
Theor. Pop. Biol. 61, 1-13. 

Guill, C, 2009. Alternative dynamical stages in stage-structured consumer populations. Theor. 
Pop. Biol. 76, 168-178. 

Hadeler, K., Van den Driessche, P., 1997. Backward bifurcation in epidemic control. Math. Biosci. 
146, 15-35. 

Hughes, T., 1994. Catastrophes, phase shifts, and large-scale degradation of a Caribbean coral reef. 
Science 265 (5178), 1547. 

Ji, Z., 2008. Hydrodynamics and Water Quality: Modeling Rivers, Lakes and Estuaries. Wiley- 
Interscience. 

Jiang, J., Shi, J., 2009. Bistability dynamics in some structured ecological models. In: Cantrell, 
R. S., Cosner, C, Ruan, S. (Eds.), Spatial Ecology. Cambridge University Press, Cambridge, 
UK, Ch. 3, pp. 181-294. 

Jordan, S., 1987. Sedimentation and remineralization associated with biodeposition by the american 
oyster crassostrea virginica (gmelin). Ph.D. thesis, University of Maryland. 

Kniskern, T., Kuehl, S., 2003. Spatial and temporal variability of seabed disturbance in the York 
River subestuary. Estuarine, Coastal and Shelf Science 58, 37-55. 

Konar, B., Estes, J., 2003. The stability of boundary regions between kelp beds and deforested 
areas. Ecology 84 (1), 174-185. 



20 



Lenihan, H. S., Micheli, F., Shelton, S., Peterson, C. H., 1999. The influence of multiple envi- 
ronmental stressors on susceptibility to parasites: an experimental determination with oysters. 
Limnol. Oceanogr. 44, 910-924. 

Liu, P., Shi, J., Wang, Y., 2007. Imperfect transcritical and pitchfork bifurcations. Jour. Func. 
Anal. 251, 573-600. 

Mann, R., Powell, E., 2007. Why oyster restoration goals in the Chesapeake Bay are not and 
probably cannot be achieved. Journal of Shellfish Research 26 (4), 905-917. 

Newell, R. I. E., 1988. Ecological changes in Chesapeake bay: are they the result of overharvesting 
the eastern oyster {crassostrea virginica)"! In: Lynch, M., Krome, E. (Eds.), Understanding the 
Estuary, Advances in Chesapeake Bay Research. Chesapeake Research Consortium, Gloucester 
Point, VA, pp. 536-546. 

Nielsen, P., 1992. Coastal Bottom Boundary Layer and Sediment Transport. World Scientific. 

Petraitis, P., Hoffman, C, 2010. Multiple stable states and relationship between thresholds in 
processes and states. Marine Ecology Progress Series 413, 189-200. 

Petraitis, P., Methratta, E., Rhile, E., Vidargas, N., Dudgeon, S., 2009. Experimental confirmation 
of multiple community states in a marine ecosystem. Oecologia 161 (1), 139-148. 

Powell, E., Kraeuter, J., Ashton-Alcox, K., 2006. How long does oyster shell last on an oyster reef? 
Estuarine, Coastal and Shelf Science 69 (3-4), 531-542. 

Powell, E. N., Klinck, J. M., Ashton-Alcox, K. A., Kraeuter, J. N., 2009. Multiple stable reference 
points in oyster populations: biological relationships for the eastern oyster (crassostrea virginica) 
in delaware bay. U.S. Fishery Bulletin 107, 109-132. 

Rietkerk, M., Van de Koppel, J., 1997. Alternate stable states and threshold effects in semi-arid 
grazing systems. Oikos 79 (1), 69-76. 

Rothschild, B., Ault, J., Goulletquer, P., Heral, M., 1994. Decline of the Chesapeake bay oyster 
population: A century of habitat destruction and overfishing. Marine Ecology Progress Series 
111, 29-39. 

Ruesink, J., Lenihan, H., Trimble, A., Heiman, K., Micheli, F., Byers, J., Kay, M., 2005. Intro- 
duction of non-native oysters: ecosystem effects and restoration implications. Annual Review of 
Ecology, Evolution, and Systematics 36 (1), 643. 

Scheffer, M., 2009. Critical Transitions in Nature and Society. Princeton University Press. 

Scheffer, M., Carpenter, S. R., Foley, J. A., Folke, C, Walkerk, B., 2001. Catastrophic shifts in 
ecosystems. Nature 413, 591-596. 

Schulte, D. M., Burke, R. P., Lipcius, R. N., 2009. Unprecedented restoration of a native oyster 
metapopulation. Science 325, 1124-1128. 

Smith, G. F., Bruce, D. C, Roach, E. B., Hansen, A., Newell, R. I. E., McManus, A. M., 2005. 
Assessment of recent habitat conditions of eastern oyster crassostrea virginica bars in mesohaline 
Chesapeake bay. North American Journal of Fisheries Management 25, 1569-1590. 



21 



U.S. Army Corps of Engineers Norfolk District, 2009. Final programmatic environmental impact 
statement for oyster restoration in Chesapeake bay including the use of a native and /or nonnative 
oyster. |http : //www . nao . usace . army .mil /DysterEIS/FINAL_PEIS/homepage . asp| 

van Rijn, L. C, 1986. Sediment transport by currents and waves. Delft Hydraulics Report H461. 



22 



